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Abstract. Reduced order models, in particular the reduced basis method, 
rely on empirically built and problem dependent basis functions that are 
constructed during an off-line stage. In the on-line stage, the precomputed 
problem-dependent solution space, that is spanned by the basis functions, can 
then be used in order to reduce the size of the computational problem. For 
complex problems, the number of basis functions required to guarantee a cer- 
tain error tolerance can become too large in order to benefit computationally 
from the model reduction. To overcome this, the present work introduces a 
framework where local approximation spaces (in parameter space) are used to 
define the reduced order approximation in order to have explicit control over 
the on-line cost. This approach also adapts the local approximation spaces to 
local anisotropic behavior in the parameter space. We present the algorithm 
and numerous numerical tests. 



1. Introduction 

The recent progresses in the numerical simulation of physical phenomena ob- 
tained through the combination of robust and accurate approximation methods 
and faster and larger computational platforms have permitted to investigate more 
and more complex problems with improved reliability. These progresses have, in 
turn, lead to new demands for the numerical simulations that are not only used 
to understand a given state but investigate control and optimization problems. In 
such problems, a large number of similar computations of a parameter-dependent 
model have to be performed associated with different values of the parameter chosen 
either randomly of guided by some well suited recursive algorithm. 

Faster solution algorithms are often not sufficient to achieve these new demands 
in many engineering applications and reduced numerical methods have been pro- 
posed as surrogates to standard numerical high fidelity approximations of given 
mathematical models. Among these methods, the reduced basis (RB) methods 
[16j [14] [15] is a class of higher-order mathematical technique that has been widely 
developed and gained in generality and reliability during the last decade. The 
greedy algorithm [18 plays a crucial ingredient and will be the main focus of the 
present paper. For parametrized time-dependent problems, the RB (greedy-based) 
strategy can be successfully combined with the proper orthogonal decomposition 
(POD) [10] as is illustrated in [8] [12]. 

The basic idea behind these reduced order methods is the notion of small Kol- 
mogorov n-width of the set of all solutions under variation of the parameter, i.e. 
that the solution manifold, embedded in a Hilbert space W, can be well approxi- 
mated with a well chosen n-dimensional subspace of W. In essence, RB methods 
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are based on a two steps strategy : the first step (off-line stage) allows to select 
particular instances of the parameters, for which a very accurate approximation of 
the solution is computed : those solutions constitute the reduced basis. In a second 
step (the on-line stage), the generic solutions (for other instances of the parameter) 
are approximated by a linear combination of these reduced basis functions. The 
interest of the approach lies in the fact that, in cases where the Kolmogorov n- 
width is small — and there are many indications for such a smallness: (i) a priori, 
such as regularity of the solutions as a function of the parameters when the set 
of parameters lies in a finite dimensional space, and (ii) a posteriori, revealed by 
a large number of simulation — these reduced model approximations require very 
few degrees of freedom and in some sense are close in spirit to other high order 
approximation methods like spectral methods, but with an improved efficiency for 
what concerns the computational complexity. 

It is well known that high order methods generally take advantage of a global 
approach by using basis functions that have a large support that, combined with 
higher regularity of the solution to be approximated (going some times up to an- 
alyticity) allows to get, with very few degrees of freedom, a very good accuracy. 
For most practical design problems in engineering though, the solution is not ana- 
lytical and most of the time regularity exists but only locally which precludes the 
interest of global approximations. This is the reason why, for instance, by breaking 
the global framework to locally piecewise global approaches, the spectral element 
method reveals superiority with respect to the plain spectral method: a trade off 
between locality and globality is generally preferred as is demonstrated in e.g. [2] for 
approximation in spatial direction by spectral element approximations. The same 
causes imply the same effects in the parameter directions for RB approximations, 
this is the reason why recently, some works have been devoted to proposing ways 
of adding locality to the reduced methods yielding parameter subdomain domain 
refinement. 

In this paper we shall focus on the certified RB framework for which, the con- 
struction of basis functions during the first stage of the algorithm is, as in most of 
the current approaches, performed through a greedy strategy based on an a poste- 
riori error estimator. For either high dimensional parameter spaces or spaces with 
large measure one may encounter the problem that the size of the reduced basis 
turns out to be larger than desired. Following the lines drawn above, first ideas in 
this context has been presented by Eftang, Patera and R0nquist [4] and further in 
[5j El [7] where the parameter space is decomposed into cells where different reduced 
basis sets are assembled. This approach presents clear advantages in the size of the 
matricial system that appears in the on-line solution procedure, and corroborates 
the natural feeling that, in order to approximate the solution at a given parameter, 
primarily those solutions in the reduced basis corresponding to parameters that 
are close to the parameter we are interested in are to be involved in the linear 
approximation. 

A drawback of the current approach [U[5j 13 [Z] however is, that in two adjacent 
parameter-subdomains, some of the parameters that are selected may be very close 
or even lie both on the common interface. This is a situation that may be en- 
countered very often when dealing with high dimensional parameter spaces, and in 
addition it is known from computations that the greedy algorithm has the tendency 
to pick sample points lying on the boundary of the parameter space. In this case, 
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the approximation on one of the subdomains does not benefit of the computation 
that was done for the other subdomain. Parameter domain decomposition may thus 
not be the ultimate approach. Another drawback of all classical greedy approaches 
is to be unable to master the size of the reduced basis and in consequence the size 
of the discrete system that will be solved in the on-line procedure. 

These reflections have motivated us to investigate the alternative discussed in 
this paper which introduces a completely new framework, which should be seen as a 
first layout of our ideas where it is unfortunately impossible to answer all interesting 
questions related to this approach without losing the big picture. 

The outline of the paper is as follows. In Section [2] we present the standard 
greedy algorithm that is widely used in reduced basis computations in order to 
have a common ground for the further discussion. Then, we introduce in Section [3] 
the concept of local approximation spaces, that is, if a reduced order approximation 
needs to be computed at a parameter point \i the approximation space is spanned 
by the N closest precomputed basis functions. Thus, the construction of the local 
approximation spaces is based on a distance function which is constructed in Section 
[4] in an empirical manner taking into account the geometry of the manifold of 
the solution set including the local anisotropics in the parameter space. Section 
[5] introduces the concept of a varying training set that is chosen in an optimal 
fashion using the possibly anisotropic distance function if the computing resources 
are limited to handle a given (and too small) number of training points. In Section 
[6] we explain some practical aspects of the on-line implementation, in particular 
how the different snapshots can be orthonormalized at the on-line stage. In Section 
[§] we present some numerical examples which uses the best approximation (L 2 - 
projection) for an explicitly given family of functions whereas an example using a 
reduced basis framework is given in Section [9j Finally, in Section [10] we draw the 
conclusions. 

2. The greedy algorithm for reduced basis approximations 

Let us first introduce a classical greedy algorithm to have a common ground to 
present our ideas. In a general setting, it is assumed that for each parameter value 
l± in a parameter domain P C R p , a /[/-dependent function v(fi) £ W of the variable 
x in a bounded domain Q C M d (e.g. d = 2 or d = 3) can be computed. The 
space W denotes some functional space of functions defined on ft and we denote 
the solution manifold by 



In order to fix the ideas this can be that v(fi) is a solution to a parameter dependent 
partial differential equation, an approximation of which can be computed by the 
e.g. classical finite element or spectral method. Let 



be the approximation space associated to the set S n. A projection operator P/v : 
W — ^ Wat is supposed to exist that can be either a parameter dependent Ritz- 
projection (reduced basis methods in case of a parameter dependent PDE) p~4j [17] 
or simply a L 2 -projection. Further, it is assumed that for each parameter value 



W P = {v(fl) £ W I /X £ P}. 



Sat • • 

be a collection of N parameters in P and 

Wat = span{i;(/x 1 ) 




■■,v(n N )} 



4 



Yvon Maday and Benjamin Stamm 



/iGPan error estimator 7y(/x;Wjv) of the approximation of v(fjb) by Pjsr(v(fjb)) can 
be computed (this can be e.g. through a posteriori analysis of the residual, see eg 
[14]). We assume that there exists two constants c > and C > such that 

cri(wW N ) < \\v(fi) - P N (v(fi))\\ <Cr/(/x;Wiv), 

for an appropriate norm on W. The case where n (/x; Wat) describes the exact error 
is thus not excluded in this general framework. 

The following sketch represent a typical greedy algorithm: Let Strain be a chosen 
finite training set Strain C P of representative points. 



ClassicalGreedy 
Input: tol. 
Output: A/", Wat, Sat. 

1. Choose (possibly randomly) fj, 1 G Strain, set Si = {jx 1 }, Wi = 
span{v(ix 1 )} and N = 1, err = max /iGHtrain 7y(/x; Wi). 

2. While err > tol 

3 . Find v N+1 = argmax /xeHtrain 77(/x; Wjv), err = max MGStrain n(w W N ). 

4. Compute v(ii N+1 ), set 8^+1 = Sjv U {m^ 1 } and set 

Wat+i = spanlWAT^^ 1 )}. 

5. Set 7V:=7V + i. 

6 . End while. 

Algorithm 1: Classical greedy algorithm. 



Remark 2.1. The introduction of the finite set S train is due to practical imple- 
mentation, because the entire scan of P is impossible: S train has to be finite, not 
too large in order the previous algorithm is not too lengthy, not too small in order 
to represent well P. Actually the definition of this finite set can evolve during the 
algorithm, see [7 . We shall elaborate on this in the framework of our approach 
later in the paper (Section^. An alternative approach is presented in PQ, where 
finding the maximum is stated as a constraint optimization problem. 

Let N denote the final size of the reduced basis space such that the accuracy 
criterion 

max 7/(/x; W^) < tol 

train 

is satisfied. Note that this final size is a consequence of the tolerance that has been 
chosen, hence, for a prescribed tolerance, we have not the hand on the complexity 
of the RB approximation method. 



3. Locally adaptive greedy approximations 

The new approach pursued in this work is in the same spirit of partitioning the 
parameter set as in [4, 5j 6, 7 . However, the difference is listed below: 

rn We do not impose a clear partition of the parameter space but rather collect 
a global set of K sample points Sk, they are preliminary selected in the off- 
line stage and the corresponding basis functions are computed and stored. 



This step is elaborated in Section 3.2 
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[a] We choose a priori the size N (< K) of the system we want to solve on-line. 
Thus, the resulting complexity of the on-line stage is user controlled. 

When a reduced basis approximation is to be computed for a certain given 
parameter value /x G Pin the on-line stage, we only use the N precomputed 
basis functions whose corresponding parameter values lie in a ball around 



/jl G P. This step is elaborated in Section |3X 
The distance function used to define the ball, containing the N closest basis 
functions, is an empirically constructed distance function that is problem 
dependent and can reflect anisotropics in parameter space and its local 
changes. The construction of the distance function is presented in Section 

m 

3.1. Local approximations spaces. As already anticipated, for a certain given 
parameter value /x G Pin the on-line stage we only use the N basis functions whose 
parameter values lie in a ball 

(l) B /1 = {AeP|d(/i,A)<r(/i)}, 

for a given distance function •)■ The radius r(/x) is tuned in such a way that 
there are actually N basis functions in the ball. Denoting by S>k the set of sample 
points for which the corresponding functions are precomputed at the off-line stage, 
the local sample space is defined by 

S /JL = B„nS K = {iieS K \iieB §A } 

with cardinality equal to Further the local reduced basis approximation space 
shall be defined by W M = span{v(/i) | jl G S M } and its associated local projection 
by P M : W — » W M in order to define the reduced basis approximation P fJL (v(fjL)) to 
v(ijl). Also in this local framework, the Galerkin setting can be used to design a 
posteriori estimates to quantify the error between v(fjb) and P M (i>(/z)), that we shall 
denote by rj(fj,] W M ). 

A similar approach is presented in [TT] where the snapshots are chosen also in the 
on-line stage but such that the approximation error is maximally reduced. Thus, 
in general, no local neighborhood is required. 

One basic ingredient of this local reduced basis approximation is thus the distance 
function the radius is adjusted to select the total number of modes N for 

each value of \i whereas the distance function can be chosen isotropic or better 
accounts for anisotropics in the parameter space. The construction of an adapted 
local distance function is a topic itself and is addressed in Section [4] Another 
ingredient is the set of sample points E>k whose construction is explained in the 
next section. 

3.2. Construction of global set of sample points. This consists of the off-line 
stage of the algorithm. The first decision is to set the number N of basis functions 
we want to get involved in the on-line procedure (i.e. the size of the matrix system 
to be solved on-line). 

A standard greedy algorithm, as described by Section [2j is then performed in a 
first stage until N + 1 basis functions are selected. The corresponding parameter 
values define the current set Sat+i and the solutions i± G §jv+i, are saved. 



■'"Note that N has not the same meaning as in the previous section any more; hence the 
cardinality of E>k is not iV but & K > N generally much larger. 
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Then, the second stage of the algorithm is performed to further enrich §jv+i- The 
construction of the further basis functions require the data of a distance function 
(as in the on-line stage). As for now we assume that the distance function is 
available, its recursive construction is explained in Section [4] (note that, by default, 
we use the Euclidean isotropic metric). We note that we can at each iteration 
and for each parameter value \i G P, compute an approximation with the N basis 
functions that lie in and estimate the error by computing 77(/L£; W M ). The next 
parameter value (s) is chosen as the one corresponding to the worse approximation 
properties. It is important to notice that using local approximation spaces, several 
basis functions can be chosen per iteration. The following procedure is proposed. 

At the end of a given iteration where we dispose of /zi, . . . , /l^x, a new iteration 
may start. The next sample point, /jlk+i, is chosen such that the estimated error 
77 (/x; Wjx) is maximized over the parameter space as in the traditional setting except 
that only a local approximation space W M is considered to compute the error (- 
estimator). Thus fiR+i = ar g max /xeH train viw Next, the maximum estimated 

error 77 (/x; W M ) is searched over all parameter values lying outside of the domain of 
influence of Hk+i- By the domain of influence of iik+i we mean 

={^P I MK+i E B^} U B^ K+1 . 

Outside of the domain of influence of Hk+i more basis functions can however be 
added. This condition assures that at most one sample point is added per ball. 
This procedure is repeated, and more and more parts of the parameter domain are 
excluded, until the remaining set is the empty set. The practical implementation 
of this procedure is given the following pseudo-code. 



Enr i chment L op 




Input : K, E>k 




Output : updated S>k 




1 . Set Strain S train , err max - — - r/(/x; W M ). 

A*t strain 




2. While Strain 7^ and err > tol 




3. Set n K+1 = argmax 77(^5 W M ), err = max 7/(/z; W„). 

A*£= Strain A*^ Strain 


4. Compute v(fi K+1 ) and set § K +i = $k U {m 




5 . Set S train := S train \{/Lfc G Strain | ^ G B^k+i or ^i K+1 




6. SetX:=X + l. 




7. End while. 




Algorithm 2: Basis enrichment loop. 





4. ANISOTROPY AND DISTANCE FUNCTION - A GENERAL FINITE DIFFERENCE 

BASED APPROACH 



In order to equip the parameter space P with an appropriate metric to Wp and 
distance function, one has to estimate, on the fly, the way the (unknown) function 
v(fi) G W depends on \i. In this manner we aim to reduce at most the number of 
elements involved in the on-line reduced basis approximation. 

Indeed, for instance assuming that the parameter domain P is a subset of R 2 , 
so that 11 = (/ii,/i2), and assuming that there exists a function ip of one variable, 
such that v(fi) = (p(^i + ^2), then the optimal selection of parameters should be 
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sought e.g. along the line [i\ = ji2 considering only a one dimensional parameter. 
Of course the fact that v has such a behavior is generally not known. This has to be 
figured out from prior computations, during the greedy algorithm. This is actually 
a quite classical quest in numerical analysis, and this valuable knowledge leads to a 
large computational reduction. As in other instances, e.g. adaptive approaches like 
the one used in finite element approximations (see e.g. [3]), the construction of such 
a distance function •) that accounts for the anisotropic behavior of changes of 
v(/j,) with respect to variations in i±. Similarly as in other approaches, the distance 
function is derived from the knowledge of the Hessian of the function v with respect 
to the parameter fi. 

We first propose a general framework to obtain an approximation of the Hes- 
sian by finite differences and then give the definition of the distance function. In 
what follows, we always assume that S train is such that its convex hull spans the 
parameter space P. 

4.1. Definition of the Hessian. The goal is to define a Hessian matrix H(/j,) for 
each point G Strain upon which the distance function will be based on. Since the 
Hessian is based on the reduced basis approximation, it is updated/constructed at 
each iteration in the off-line greedy algorithm (Section |3.2| ). To do so, compute a 
reduced basis approximation PuM aO) °f V (/ Jj ) based on the current local reduced 
basis, as explained in Section [37l[ on the stencil /x(a) = /jl + Sf=i c^/A w ^h 
M £ Strain and where 8jA l , i = is a positive small increment in the i th 

direction in the parameter and a 1 = —1,0, 1, such that, at most, two of them are 
non zero. Let {y>n}n=i denotes the basis of W M used to build P /Lt (^(/x)) such that 

N 



From the 3 P approximations P M (i>(/x)) of v(fi(a)) on the surrounding stencil /x(a) 
of /j,, one can, by finite differences, define a Hessian matrix with value in R as 
follows 

N 
n=l 

where 

v n (fi + Sfi z + Sfi J ) - v n (fjL - Sfi z + Sfi J ) - v n (/jb + Sfi z - 5{i J ) +v n (fi- 5{i z - 5{i J ) 

= Wsjp • 

The term DijV n (fi) describes the anisotropy of the n-th mode while v n (fi) indicates 
its weight in the approximation. 

Remark 4.1. In the case one is only interested in a scalar output functional s : 
W — >• R of the solution, i.e. s(v(/i)) approximated by s(i? rb (/x)) ; one can also 
consider the Hessian do be defined as 

Hij(fjb) = Dijs(v Ih (fjb)). 

Remark 4.2. The current construction of the Hessian based on a finite difference 
scheme involves quite a large number of points when the dimension p becomes too 
large. Alternative exist based only on the only points in S>k that need still to be 
investigated and improved. 
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At this level we dispose for any /x G P of a method to compute Hij(fji). It 
remains to explain how we construct the metric tensor M from the knowledge of 
the Hessian matrix H. Indeed, given \i G P, we proceed in a standard manner and 
perform an eigenvalue decomposition 

H(fjb) = V AV T 

where V is an orthonormal matrix and (A)^ = A$ a diagonal matrix consisting of 
the eigenvalues A;. Consider the diagonal matrix 

\Mu = M> i = >P> 
and the associated symmetric semi-positive definite matrix M(/x) = V |A| V T to 
define the metric tensor. 

4.2. Construction of the distance function. The construction of the distance 
between two parameter points /xi, /12 G P, denoted by /12), in the metric space 
given by the metric tensor M is ideally defined by the length of the geodesic curve 
linking fii and M2- Indeed, if 7 : [0,1] — >> P is a parametrization of the geodesic 
from i±\ to i±2 we would like to compute 

jT 1 yJh'(t)) T M{-r{t))y{t)dt. 

The above integral is however not computable due to the lack of knowledge of 
the geodesic curve. The geodesic curve itself could be approximated, but this 
would slow down the on-line computation drastically. We therefore propose to 
approximate the above integral roughly by the following trapezoidal rule to define 
the distance 

d(Mi>M2) = ~ Mi) T O2 - Mi) + \\J(Vi ~ Vi) T M (^) (M2 - Mi) 

where we also took the approximated 7' '(£) w (//2 — Mi) for t = 0, 1 (which is correct 
when dealing with the Euclidian metric). 

We note that the construction of the Hessian H(fjb) and thus the distance function 
depends on the approximation P^v and thus on Sk- That is, for each reduced 
basis model we can define a distance function and we recursively update the distance 
function at each step of the stage 2 of the greedy algorithm. Therefore, the distance 
function is, simultaneously with P fJL (v(fjL)), converging to a final distance function 
that is used in the on-line procedure. 

Let us make the following comments: 

H~| It is easy to see that our distance function is symmetric by construction. 
The matrix M(/x) is semi-definite positive. Allowing for zero eigenvalues 
is indeed a desired property since it allows to shrink certain directions in 
parameter space where there is no variation in the solution with respect to 
the parameter. 

1~3~1 Defining a mesh consisting of p-dimensional simplices with vertices being 
the set S train allows to interpolate/represent quantities only defined on 
Strain by means of piecewise linear functions on the mesh for any parameter 
value /1 G P, since the convex hull of S train coincides with the parameter 
space P. We can therefore easily access the interpolated metric tensor 
M(/x), for any fi G P, by means of its interpolated component functions 
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Remark 4.3. As an alternative to interpolating the coefficients in a p- dimensional 
space one can construct the Hessian on-line as explain above for any fj, G P. 

5. Adaptive training sets based on Hessian 

The goal of this section is to present an approach that further allows to reduce 
the off-line computational cost. Here, we focus on the training set S train . It is 
aimed to keep its cardinality as small as possible, but large enough to capture 
the local geometry of the parametrized system. We aim to give an answer to the 
following question: Given the computational resources of handling Qm points in 
Strain, how do we place them optimally in the parameter space? 

We propose to place the Qm points uniformly in P with respect to the empirically 
constructed metric tensor M(/x) /r(fj,) 2 . The factor l/r(/L£) 2 is motivated by the fact 
that a ball defined in ([T| corresponds (approximatively) to a unit ball in the metric 
defined by M(/x)/r(/zp. Since the metric gets updated at each iteration of the 
greedy algorithm, we also need to adapt the training points S train at each iteration. 
Further, we propose that the cardinality of S train is an increasing function of the 
inverse of the actual error (-estimation) err. Let Q m be the minimal cardinality 
of S train at the beginning of the greedy algorithm and Qm the maximal target 
cardinality at the end of the algorithm once the tolerance criteria err < tol is 
reached, then we define 



Q(err) 



log(err) + Q r , 



log(tol) 

to be the number of points at the next iteration in S train . 

To summarize, the framework that we propose is that the training set Strain of 
cardinality Q(err) is constructed such that the points are uniformly distributed in 
the empirically constructed metric tensor M(/x) /r(/j,) 2 . 

Our implementation of this concept is based on the Delaunay mesh defined by 
the points of S train . We use the software FreeFem++ [9] which provides a uniform 
mesh with respect to M(/x). Other implementations are of course also possible. 

Remark 5.1. Enriching the training set during the sampling process is not a new 
idea. An adaptively enriching greedy sampling strategy was proposed in [7] however 
it is not based, as it is the case in the present paper, on a learning process of the 
geometry of the manifold. In addition our approach is more precise since it not 
only detects where the manifold of the solutions (as a function of the parameters) 
is complex but also the directions that are important to sample and those that are 
not. Second in our approach the size of the discrete system to be solved on-line is 
determined by the user leading to a controlled simulation cost. 

6. On-line orthonormalization of basis functions 

In the standard approach, it is well known that the stable implementation of 
the reduced basis method requires the orthonormalization of the snapshots 
\± G §tv [13]. This orthonormalization is prepared off-line and the vectorial space 
Wat spanned by the series \v[\i) \ \i G Sat} is not affected by this change of basis. 

In our alternative method, there are a large number of approximation spaces, 
almost one for each /j, : W M = span{v(/i) \ p, G S M }. We cannot compute all various 
orthonormalized basis sets for all possible W M off-line. This would be too costly 
and using way too much memory. 
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Nevertheless, we can still prepare the orthonormalization process so that, the 
on-line orthonormalization is feasible. For any f±i and fij in let us compute 
off-line the scalar products 

(2) Mij = (vfa^vik)), i,j = l,...,K 

where (•,•) represents an appropriate scalar product : say either L 2 or H 1 type. 
Note that we can even compute only those scalar products corresponding to pairs 
of parameters (/x^, fij) that are close, indeed if fti and fij are distant in P, they will 
never be in the same S M for any /jl G P. 

In the on-line stage, once the approximation in W M will be required, we can 
orthonormalize the basis set {v(fi) | fx G S M } first by declaring an order in the 
parameters 

(3) {fx G S M } = {mi, ^2, Mat} 

then perform a classical Gram-Schmidt orthonormalization process 

Ci = /3iiu(/zi), 

( 2 = P22V(V2) + /?2lCl, 
iV-1 

(4) Civ = Pnnv(^n) + 

2=1 

where the coefficients /3 n i are chosen so that ( n is orthogonal to 

span{v(/L^) | i = 1, .., n — 1} = span{^ | i = 1, .., n — 1} 

and the coefficients /3 nn are chosen so that £ n is norm 1. It is well known that 
these coefficients exist, are unique and that, in order to compute them, we have to 
know each scalar product (v(fj, n )Xi), i = 1, — 1 that can be easily computed 
mimicking the Gram-Schmidt procedure Q but using the precomputed coefficients 
Mij. For the construction of the coefficients /3 n j, we proceed as follows (it helps 
to have in mind that pkn stands for — (v(fJLk),Cn)) : for all n = 1, . . . , TV apply 
recursively 

n-1 

nn — y ] \Pnj | 
J = l 

Pnn — i 

a 

finj ' ~ •> k = 1. . . . . n — 1, 

a 

n-1 

(5) hn = -PnnM nk + ^2P nj j3 kj , fe=l,...,JV. 

Note that the construction of a is not subject to instabilities since the sum involves 
only positive values. 

Next, we notice that the basis functions {( n }n=i can alternatively be expressed 
by the following linear transformation 

Ci = 7ii^Oi), 

C2 = 722^(^2) + 721^1), 
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N 

(6) Cn = ^2~f N iv(in), 

2 = 1 

from the set of basis functions {v(fjL n )}^ =1 . It is our final intention to construct 
the coefficients j n i representing the change of basis functions from {v(fjL n )}^ =1 to 
{(n}n=i as defined by ([6| in a stable way and without any operation depending on 
the length of the vectors representing v(fi n ) to guarantee stability and good on-line 
performance. 

Remark 6.1. Note that the coefficients 7^ could be obtained based on (|6| directly. 
Indeed, each ( n is orthogonal to span{v(ni) \i = 1, .., n — 1} allowing to compute the 
coefficients j n i, i = 1, . . . ,n — 1. This however involves solving a linear systems 
based on the mass matrix with coefficients = (v(fjbj), v(fJLi)), i,j = l,..,n — 1, 
which can be ill conditioned due to the fact that the set of functions {v(/Jti)}i is 
heavily linearly dependent. 

The set of coefficients {j n i}ni=i can be obtained by 0(N 3 ) operations from the 
coefficients /3 n j, through a triangular process. Indeed, the coefficients 7^ can be 
obtained by the following recursive formula 

{YJk=i Pnklki if i < n, 
p nn if i = n, for all n = 1, . . . , N. 

if i > n, 

As mentioned above, the set of basis functions {Cn}n=i l ea ds to well-conditioned 
matrices, the precomputed quantities are nevertheless expressed in the basis {v(n n )}, 
but the change of basis functions is described by the coefficients {7m}^=i as stated 
in equation ([6|. 

As an illustration assume that the reduced basis problem to be solved is about 
the approximation of the solution to a parametrized PDE, written in a variational 
form as : for a given parameter value /u, find G W M such that 

a(u(fjb), v; 11) = f(v; /x), Wv e W M , 

where it is assumed, for the sake of simplicity that a and / are affine decomposable 

Qa Qf 

(8) a(w,v,fi) = ^2g q (fji)a q (w,v), f(v^)^h q (^)f q (v). 

q=l q=l 

The ultimate goal is to construct the well-conditioned stiffness matrix a(Q , l*) 
derived from the off-line pre-computation of the series A\- = a q (v(/j,j),v(/j,i)), for 
1 < h j < N. To do so, follow the next steps: 

H~| Compute the coefficients {/?m}rM=i following ([5|. 

Compute the coefficients {7m}^i=i following (|7|). 

[3] For each q = 1, . . . , Q a , the matrix a q (Q, d) can De computed by a q (Q, Q) — 

Build the solution matrix by summing a(QXi'-> m) = J2 q =i dqif 1 ) a q((ji(i)- 
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Similarly, the evaluation of each component of the vector /x), 1 < i < N can 
be performed in 0(QfN) operations, and the inversion of the associated matricial 
problem is done in G(N 3 ) further operations. 

7. Overview of the framework 

Before we start with the numerical tests, an abstract description of both the 
off-line and on-line parts of the locally adaptive greedy algorithm is given in the 
following box: 




========= Stage 1 ========= 

1 . Perform a classical greedy algorithm to select N + 1 basis functions. 
========= Stage 2 ========= 

2. For each fi G S train , do (parameter space search) 

a) Compute the error estimate r/(/x,W M ). 

b) Update the Hessian matrix Hij(fjL). 



4. Enrich the set of basis functions (Enrichment Loop, Section 3.2) 

5. Create a new training set Strain (Section [5|. 

6. Go to 2. until tolerance tol is achieved. 

Algorithm 3: Locally adaptive greedy algorithm (Off-line). 



LocallyAdaptiveGreedy On-line 
Input: K, Wk, \i 
Output: P M M//)) 



Find the approximation space 

a) Compute the distance <i(/x, jl) for each jl G Sk (Section 4.1). 

b) Choose the N closest jl G §k to form § M and ^ 
span{v(/i) | fx G § M } (Section 3.1). 



2. Orthonormalize the basis functions \v{fi) \ fx G S M } to build a new basis 
{£i 5 --- 5 £iv} (Section [6). 

3. Compute the projection P^{v(fi)) based on the basis . . . ,£at}. 

Algorithm 4: Locally adaptive greedy algorithm (On-line). 



8. Numerical results using the L 2 -projection 

Before we pass to a full Reduced Basis model based on a parametrized Galerkin- 
projection, we use a parametrized family of explicitly given functions combined 
with the L 2 -projection to define the projection operator Pjy in order to test our 
framework. Reminding Cea's lemma which is a consequence of the Galerkin ap- 
proach, studying the best-approximation is sufficient to test the method of locally 
adaptive approximation spaces. The error is computed in an exact manner and no 
a posteriori estimation is used. 

Further, in a first step, Section [8TT[ we only test the feature of local anisotropic 
approximation balls without adapting the training set using the distance function. 
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That is we use the same fixed Strain throughout all steps of the algorithm. In a 
second step, Section [8^2} we test the entire algorithm. 

For all numerical tests, Q and P are discretized by a regular lattice of 75 x 75 
points and the target tolerance is set to 10 -4 if not otherwise mentioned. 

8.1. Numerical results with fixed training set. 

Test 1. We start with presenting a numerical example to illustrate the benefit of 
the local anisotropic approximation spaces. Consider the function 

, , f x _ ^ [ -(xi -Q.l(/ii -/i 2 )) 2 _ (x 2 - (/ii +/i 2 )) 2 " 

x e n = (-1, l) 2 , /z E P = [-0.5, 0.5] 2 

that exhibits a constant anisotropy of parameters over the whole parameter space. 
The dependency in the (/ii+/i2) -direct ion is ten times stronger than in the (fii—^)- 
direction. Figure [I] (left) illustrates the evolution of the number of basis functions 
that are necessary to be computed versus the achieved accuracy in the L°°-norm 
for N = 20 for the anisotropic approach compared to the same algorithm but using 
isotropic local approximation spaces as comparison. One can observe that using 
the anisotropic approach is beneficial in terms of the number of truth solutions to 
be computed at the off-line stage. Figure [l] (right) presents the same quantity but 
for varying N = 20, 30, 40. Not surprisingly the number of needed basis functions 
decreases while increasing N and more and more the exponential convergence estab- 
lishes as is manifested for the standard greedy algorithm in this example. We can 
also observe that for a lower N more basis functions can be included per iteration. 

Figure [2] illustrates the local approximation spaces expressed in P, the radius 
as a function of the parameter and the sample set S>k for the anisotropic and the 
isotropic approach. 




Figure 1. Test 1: Accuracy with respect to the number of truth 
solutions to be computed in comparison with the isotropic ap- 
proach for TV = 20 (left) and the for different values of N = 
20,30,40 (right). 
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Local approximation spaces 



Sample points. 



Figure 2. Test 1: Local approximation spaces for selected param- 
eter values (left), radius as a function of the parameters (middle) 
and sample points (right) for anisotropic approach (top) and the 
isotropic version (bottom, as comparison) with TV = 20. 



Test 2. The previous example is in some sense idealized since the anisotropy is not 
varying throughout the parameter domain. In this regard, the next example 



/2 («;//) = exp 



(xx-(/i?+M2)) 2 (^2-(/i?+^2)) 2 



0.01 



0.01 

x e n = (-1, l) 2 , ^ e P = [-0.5, 0.5] 2 

is more interesting. It presents a family of parametrized functions where the func- 
tions (as functions of x) are constant along concentric circles around the origin in 
parameter space. As an example, all four corners of the parameter space define the 
same function. 

We compare again the number of required truth solutions for the anisotropic 
approach compared with the presented algorithm but using isotropic approximation 
spaces, which is presented in Figure H (left) for TV = 5. We observe that 60% of the 
number of truth solutions to be computed can be saved by constructing the problem- 
dependent approximation spaces. The corresponding local approximation spaces, 
radii and sample points S are illustrated in Figure [4j Figure [3] (right) also presents 
the number of required truth solutions for N = 5 but for different tolerance levels. 
We observe that starting from a critical value N c of computed truth solutions, that 
depends on the tolerance, the curves show different behavior. During the iterations 
before N c the complete parameter region has to be scanned because the tolerance 
is achieved nowhere in the parameter space. During the iterations after 7V C however 
the tolerance is achieved in an increasing region of the parameter space and this 
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part needs no longer to be not scanned. The enrichment of new truth solutions is 
then limited in space and the final phase of the greedy algorithm converges faster. 

Additionally we present in Figure [5] the local approximation spaces for the dif- 
ferent tolerance levels. One can clearly see how the approximation spaces take into 
account the local geometry that also allow non-convex approximation spaces. 




J I , I I I ! I I I I , I I I ! I . I F I , , , | I | , , | I | , | | I | , | | I , , , , I , 

200 400 600 800 100 200 300 400 500 

Number of truth solutions computed Number of truth solutions computed 



Anisotropy vs. Isotropy Different tolerances 

Figure 3. Test 2: Accuracy with respect to the number of truth 
solutions to be computed in comparison with the isotropic ap- 
proach (left) and the for different end tolerances (right). In both 
cases there holds N = 5. 



Test 3. The next example 

(x 1 - (iii + 3/i 2 )) 2 (x 2 ~ (3/ii - /i 2 )) 2 



f 3 (x-fi) = exp 



0.1 + 5|/ii + 3/i 2 | 



0.1 



tie 



-0.5, 0.5] 2 



5|3/ii - /i 2 | 

x g n = (-i,i) 2 

is interesting in the sense that it presents an almost singularity in parameter space at 
the origin. The performance of the anisotropic approach, compared to the isotropic 
version, is presented in Figure [6] Again, the anisotropic approaches outperform the 
isotropic one and about a 56% of computations of truth solutions can be saved. 

The local approximation spaces, the radius and the sample points with TV = 10 
are presented in Figure [7j We observe a posteriori that the training space S train was 
not sampled fine enough. Indeed in the region around the origin (and the cross for 
the isotropic version), every training point is included in the set of sample points. 
This also explains the sudden drop of the convergence in Figure [6| in particular for 
the isotropic version. This is a known difficulty in greedy methods and our approach 
with adaptive training sets presents a solution to this problem. The corresponding 
numerical results are presented in the next section. 

Further, we plot the local approximation spaces for different values of the toler- 
ance levels in Figure [8j We observe that the originally non-connected approximation 
spaces are becoming connected with increased tolerance. 

We also recognize that the scheme detects the cross where the behavior of the 
function is most singular and present in Figure [9] the sample points that were chosen 
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Figure 4. Test 2: Local approximation spaces for selected param- 
eter values (left), radius as a function of the parameters (middle) 
and sample points (right) for the presented approach (top) in com- 
parison wit the isotropic version (bottom) for N = 5. 





0, 


f 


0.3 






i 


••• / 


-0, 


{ • ) 

V / 



^1 ^1 ^1 

Tol=10- 3 Tol=10- 4 Tol=10- 5 



Figure 5. Test 2: Local approximation spaces for different toler- 
ance levels. 



for different functions of the type 



h,tu&( x ; /•*) = ex P 



0.1 + 5|6| 



0.1 

x ett = 



f 5|6| 
(-1,1) 2 



-0.5,0.5] s 
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where £1 and £2 are functions of = (/ii,/i2) using an error tolerance of 10 -3 in 
this case. 




1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 

200 400 600 800 

Number of truth solutions computed 



Figure 6. Test 3: Accuracy with respect to the number of truth 
solutions to be computed in comparison with the isotropic ap- 
proach for N = 10. 




Local approximation spaces 



Radi 



Sample points 



Figure 7. Test 3: Local approximation spaces for selected param- 
eter values (left), radius as a function of the parameters (middle) 
and sample points (right) for the presented approach (top) in com- 
parison wit the isotropic version (bottom) for N = 10. 




Figure 8. Test 3: Local approximation spaces for different toler- 
ance levels. 




£l — M1+3/X2, £2 = -M2 £l = Ml + M2, £2 = Ml - M2 

Figure 9. Sample set patterns for different type of parametrized 
functions with singular behavior on two crossing line with different 
slopes. 



8.2. Numerical results with adapted training sets. In order to compare the 
accuracy of this version with adapted training sets we introduce a fixed test set of 
sample points S test . In the following examples it consists of a lattice of 75 x 75 
points that corresponds to the fixed training set of the tests of the previous section. 
We use the command adaptmesh in FreeFem [9J in order to adapt the sample points 
by constructing a desired uniform mesh in the new metric defined by the Hessian. 

The target advantage of our new approach with an adapting training set is not to 
use fewer computed solutions than the approach using a fixed training set, at least 
in low parameter's dimension but rather provide an improved reliability. Indeed 
a non-adapted training set may leave large errors in parameter's regions that are 
not sufficiently sampled by the fixed training set. On the contrary the adapted 
set tracks these unexplored parameter's regions where the RB approximation error 
remains large. The points of the training set are chosen wisely using the information 
available from the local geometry of the system. Note however that for problem 
where the manifold Wp happens to be uniformly regular, the benefit of our approach 
will be on the computational efficiency as a result of the use of a small training set 
in the beginning that is gradually increasing. 
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Figure [10] illustrates the total number of basis functions that are necessary to be 
computed versus the achieved accuracy of the proposed algorithm with adaptive 
training sets for the test cases 1, 2 and 3. In each plot, we present the error of the 
version using an adaptive training set on the adaptive training set and the error 
evaluated on the fixed testing set S test . The accuracy of the version with a fixed 
training set (which equals S tes t) is given as comparison. 



Test 1. From Figure 10 (left), one can observe that the error on the training set 
and on the test set of the approach with an adaptive training set are similar. This 
illustrates that the accuracy is not only guaranteed on the adaptive training set 
but also satisfied on the test grid. Both errors are also similar to the error of the 



algorithm using a fixed training set as presented in Section 8.1 

The adaptive version is as good as the version on a fixed training set, but uses 
less evaluation of the error estimate as is shown in the following table: 

Number of training point evaluations using adaptive training sets: 150' 923 
Number of training point evaluations using fixed training sets: 455'625 

We observe a gain of 65% fewer error evaluations. 



Figure 11 presents the training set, local approximation spaces, radii and sample 



points for this calculation. 




Test 1: TV = 20 



Test 2: N 



Test 3: TV = 10 



Figure 10. Accuracy with respect to the number of truth solu- 
tions that need to be computed for all three test case. 



Test 2. Next, we consider the second example of the previous section. Figure [10 



(middle) compares this version (with adaptive training sets) with the version of 
the previous section using a fixed training set. One can observe again that the 
performance is similar than the version using a fixed training set, and that the 
accuracy is also satisfied on the fixed test set S test . 

Comparing again with the version using a fixed training set, the version with an 
adaptive training set requires less evaluation of the error estimate as is shown in 
the following table: 



Number of training point evaluations using adaptive training sets: 60'202 
Number of training point evaluations using fixed training sets: 151 '875 
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Training points. Local approximation Radius. Sample points, 

spaces. 



Figure 11. Test 1: Training points, local approximation spaces 
for some specific parameter values, radii and sample points for 
greedy algorithm with adaptive training set and N — 20. 



We observe a gain of 70% fewer error evaluations. 



Figure 12 presents the training set, local approximation spaces, radii and sample 
points of this calculation. 




Training points. 




Local approximation 
spaces. 





Radius. 



Sample points. 



Figure 12. Test 2: Training points, local approximation spaces 
for some specific parameter values, radii and sample points for 
greedy algorithm with adaptive training set and N = 5. 



Test 3. Finally, we apply the new algorithm with an adaptive training set to the 
third example. As already mentioned, the solution of this example shows very 
strong changes with respect to the parameter and the version with a fixed training 
set does not resolve this since the training set was not sampled fine enough. We 
expect the version with the adaptive training set to perform accurately on the 
entire parameter space, and thus in contrast using more sample points than the 
version with a fixed training set. Figure 10 (right) presents the evolution of the 
number of sample points versus the accuracy of the algorithm. One can observe 
that this approach (with adaptive training set) uses indeed more sample points than 
the version on a fixed training set as explained above. Finally, Figure 13 presents 
the training set, local approximation spaces, radii and sample points. We observe 
that the training points are in accordance with the sample points, and that a more 
dense sampling is indeed required around the origin. We observe that the local 
approximation spaces are now connected also for the tolerance of 10 -4 . 
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In order to illustrate the benefit of using adaptive training sets also in this case 
we consider a test sample of 100 x 100 uniformly distributed points in the region 



0.05, 0.05] 2 around the origin. Figure 14 illustrates the error distribution using 



the online procedure generated using a fixed (left) and an adaptive (right) training 
set. The maximum error is 0.046 resp. 1.29 • 10 -4 . While the error tolerance is 
almost satisfied in the latter case, it is clearly not the case for the former approach. 
Further, the number of error evaluations is given in the following table: 



Number of training point evaluations using adaptive training sets: 181 '672 
Number of training point evaluations using fixed training sets: 247'500 



Thus, the adaptive version uses still less error estimator evaluations and is more 
accurate in the region around the origin. Also, the error is more equally distributed. 




Training points. Local approximation Radius. Sample points, 

spaces. 



Figure 13. Test 3: Training points, local approximation spaces 
for some specific parameter values, radii and sample points for 
greedy algorithm with adaptive training set and N = 10. 





^ 1 ^ 1 



Fixed training set. Adaptive training set. 

Figure 14. Test 3: Error distribution on a region around the 
origin of the greedy version using a fixed (left) and an adaptive 
(right) training set. 
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9. Numerical results using a Galerkin-projection 

Test 4- Finally, we present an example involving the full reduced basis framework. 
We consider a two-dimensional steady convection-diffusion equation on Q = (0, l) 2 . 
Homogeneous Dirichlet boundary conditions are imposed on the left, upper and 
right boundary, i.e. on To = {(x,y) G dft \ x = 0, or x = 1, or y = 1}. On the 
lower boundary, i.e. on T g = {(x,0) G dVt \ x G [0,1]}, we impose the following 
Dirichlet data 

u(x, 0) = g(x) := x X[o,o.5) + (-40 x + 21) X[o.5,o.55] + - 1) X[o.525,i] , 

(essentially a slightly regularized saw tooth function) where xx denotes the charac- 
teristic function for the mono-dimensional region X. Thus, the boundary condition 
presents a strong gradient between (0.5,0) and (0.525,0). Next, we introduce two 
parameters, the diffusion coefficient e = 10 Ml and the direction of the convection 
parametrized by an angle /i 2 : 

f3(fi 2 ) = (sin(/i 2 ),cos(/i 2 )) T . 

Denoting the parameter vector by /j, = (/ii,/i 2 ) G P = [—2.5,0] x [— 7r/4, 7r/4], the 
parametrized problem writes as: for any /x G P, find u(fjb) G if 1 (ft) such that 

(9) -lO^ 1 Au(/x) + V • (/3(/x 2 )w(/x)) = 0, in ft, 

(10) u(fj,)=g, on T 9 , 

(11) u(/i) = 0, onr . 

The solutions presents a convected near discontinuity, that is also subject to a 
diffusion process, along a parametrized velocity field and parametrized diffusion 
coefficient. To solve the truth problem, we use a finite element solver using a 
uniform quadrilateral mesh and tensorized polynomials of degree 3 with mesh size 
h = 0.025. Different solutions corresponding to different parameter values are 



illustrated in Figure (15) and we highlight that the discretization is fine enough in 
order to avoid strong oscillations due to under-resolution related to the dominant 
convection for all parameter values. For the error computation/estimation that is 
used in the greedy algorithm, we use for sake of simplicity the exact error between 
the truth solution and the reduced basis approximation. 

This problem can be cast as a variational problem with bilinear- and linear forms 
that are affine decomposable as described in ^ with 

ai{w, v) = / VwVvdx, 9i{l*) = 10 Ml , 

a 2 (w,v)= / d x wvdx, 02 (a0 = sin(/i 2 ), 

a 3 (w,v)= / dywvdx, g 3 (fjb) = cos(/i 2 ). 

The force field is equal to zero in this case, however in the practical implementa- 
tion we impose the boundary condition nodally and project the equation onto the 
nullspace of the (discrete) Dirichlet trace operator. In this manner, the boundary 
condition appears as a right hand side. Indeed, the weak form of ([9| can be stated 
as: for any \i G P, find u(fi) G Hq(Q) such that 
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where u(fji) = u(fJb) + u g and u g \r g = g. Thus we have a modified linear form 
f(v]fjb) with an affine decomposition inherited of the one from the bilinear form, 
i.e. 

3 

f(v-fi) = -a(u g ,v-fi) = - ^2g q (fi) a q (u g ,v). 

q=l 



Figure 16 illustrates the local approximation spaces for some specific parameter 
values, the radius as a function of the parameter and the chosen sample points for 
N = 20 for the version with a fixed training set. We can observe that the symmetry 
of the problem is reflected in all three quantities and that the shape of the local 
approximation spaces does vary. Further, in Figure [16] we illustrate the training 
set, the local approximation spaces for some specific parameter values, the radius 
as a function of the parameter and the chosen sample points for N = 20 for the 
version with an adaptive training set. 

In Figure [l8| (left), we present the evolution of the error with respect to the num- 
ber of computed truth solutions during the locally adaptive greedy algorithm for 
N = 15, 20 and also plot the convergence of the standard greedy algorithm, which 
of course converges faster but with a total number of basis functions equal to 121 at 
the on-line stage. In Figure [l8| (right), the evolution of the accuracy on the training 
set of the version with the fixed and the adaptive training sets are compared. The 
value of about 4- 10 -4 of the L°°(P)-norm on the test grid is surprisingly inaccurate. 
Further numerical investigations have shown that the violation of the tolerance cri- 
terion is mainly focussed on the lower and upper boundary of the parameter space. 
The problem on the lower and upper boundary is caused how our implementation 
of generating the adaptive training set using FreeFem++ is placing the boundary 
points. Indeed, the criterion is violated on only 0.35% of all points of the test set. 
If one studies the error on a test grid of 87 x 99 points on [—2.5, 0] x [—0.7, 0.7] the 
error of the version with a fixed training set performs with an error of 1.13 • 10 -4 
and the version with an adaptive training set shows an error of 1.06 • 10 -4 . 

We finally present in Table |??| some computation times of the on-line stage and 
compare it to a standard greedy-approach in order to highlight possible improve- 
ments that have to be made in a future work. As can be seen in Algorithm 4, the 
on-line stage consists of three steps: 1) Finding the N closest sample points, 2) 
effecting the on-line ortho-normalization and 3) solving the linear system. In Table 
the average computing times over a sample of lO'OOO random parameter values 



are given and the computations are done in Matlab restricting the computations 
to a single thread using a MacBook Pro with an Intel Core i7 processor with 2.66 
GHz. We observe that the proposed locally adaptive greedy algorithm for N = 20 
basis functions is significantly slower by about one order of magnitude than the 
standard greedy algorithm using N = 121 basis functions. While solving the linear 
system is of course faster for 20 unknowns than for 121, most time is spent by 
finding the closest sample points and by effecting the ortho-normalization. Certain 
improvements can certainly be done respecting the implementation of the proposed 
algorithm such as building a tree-like structure in the parameter space that indi- 
cates the considered sample points that might be used in each box. This tree can 
be constructed off-line. At the on-line stage and for a given parameter value /j, one 
has to identify the box that contains /j, and then searching for the N closest basis 
functions only among the sample points proposed by the selected box. Further, the 
triangular nature of the matrix 7 in the on-line ortho-normalization process can be 
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exploited. In addition, the error estimation/certification is the most expensive part 
in most cases of the on-line stage of the reduced basis method. These computation 
times have not be accounted for since we worked with the exact error computation. 
Summarizing, these computation times indicate possible bottlenecks of the present 
approach but do not show the entire truth. 




\x = (-2.5, -tt/4). /x = (-2, 0). /x = (-1, tt/4). 



Figure 15. Test 4: Truth solutions for different parameter values 
which consist of the diffusion coefficient and the direction of the 
convection field. 




Local approximation spaces. Radius. Sample points. 



Figure 16. Test 4: Local approximation spaces for selected pa- 
rameter values (left), radius as a function of the parameters (mid- 
dle) and sample points (right) for N = 20 for the version using a 
fixed training set. 



10. Conclusions 

We presented a framework that can be seen as a generalization of the classical 
greedy algorithm that is widely used in reduced basis methods. The presented al- 
gorithm introduces local approximation spaces (in the parameter space) that also 
account for local anisotropic behavior instead of a global approach. The key idea 
was to consider the N closest basis functions (for a fixed N) where the distance 
is measured in an empirically built distance function which can be constructed on 
the fly and which is problem-dependent. In consequence, the number N of basis 
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Training points. 




Local approximation 
spaces. 



Radius. 



Sample points. 



Figure 17. Test 4: Adaptive training set, local approximation 
spaces for selected parameter values, radius as a function of the 
parameters and sample points for N = 20 using the version using 
an adaptive training set. 




....... N = i 5 

---- N = 2 o 

— ■■ — Standard greedy algorithm 



Number of truth solutions computed 
Using a fixed training set. 



- Error on training set (with adaptive training set) 
Error on test set (with adaptive training set) 
Greedy on fixed training set (as comparison) 




Number of truth solutions computed 
Using an adaptive training set. 



Figure 18. Test 4: Accuracy with respect to the number of com- 
puted truth solutions for the version with fixed training set (left) 
for TV = 15, 20 and the version with an adaptive training set (right) 
for N = 20. 



functions considered for an approximation is user controlled. Further, it is pro- 
posed to optionally place the training set uniformly distributed in this empirically 
constructed (approximative) metric space. Numerical tests illustrate the charac- 
teristics of this approach. In our future work, the emphasis will be shed on the 
implementation of this algorithm for high dimensional parameter spaces, the main 
challenge in this direction will be to design properly a set of equi-distributed points 
with an adapted non isotropic metric (in particular our current implementation). 
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[ms] 


Total on-line 
time 


Search of TV- 
closest sample 
points 


Ortho- 
normalization 


Solving the lin- 
ear system 


Fixed training 
set 


2.9375 


2.244 


0.59266 


0.10083 


Adaptive train- 
ing set 


2.9182 


2.2285 


0.58874 


0.10091 


Standard 
greedy algo- 
rithm 


0.3977 








0.3977 



Table 1. Test 4: Average on-line computation times over a ran- 
dom set of lO'OOO parameter values. 



IN HIGH DIMENSIONS" held in June 2011, organized with the support of the 
Fondation de Sciences Mathematiques de Paris. 
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